function G = gabor_fn2(a, S, K, sigma, lambda, psi, gamma)
%GABOR_FN2 Matlab code from Wikipedia.
% a      : a > 1
% S      : number of scales.
% K      : number of orientation, theta = n/K * pi.
% sigma  : sigma of Gaussian envelope
% gamma  : spatial aspect ratio.(Shape of Gaussian envelope)
% lambda : wavelength of the sinusoidal factor.(1/f0?)
% psi    : pharse offset

% theta  : orientation of the normal to the parallel strips of a Gabor
%          function

sigma_x = sigma;
%sigma_x = r/f;
sigma_y = sigma/gamma;
f = 1/lambda;
nstds = 4;
%a = (Uh/Ul)^(-1/(S-1));
sigma_u = 1/2*pi*sigma_x;
sigma_v = 1/2*pi*sigma_y;
%sigma_u = (a-1)*Uh/((a+1)*(2*log(2))^0.5);
%sigma_v = tan(pi/(2*k))*(Uh-2*log(2)*sigma_u^2/Uh)*(2*log(2)-(2*log(2))^2*sigma_u^2/Uh^2)^(-0.5);
c0 = 1/(2*pi*sigma_x*sigma_y);
c1 = 2i*pi*f;
G=cell(S,K);
for n = 1:K
    theta = n * pi / K;
    % Bounding box
    xmax = max(abs(nstds*sigma_x*cos(theta)),abs(nstds*sigma_y*sin(theta)));
    xmax = ceil(max(1,xmax));
    ymax = max(abs(nstds*sigma_x*sin(theta)),abs(nstds*sigma_y*cos(theta)));
    ymax = ceil(max(1,ymax));
    xmin = -xmax; ymin = -ymax;
    [x,y] = meshgrid(xmin:xmax,ymin:ymax);
    % Rotation
    x_theta_r = x*cos(theta)+y*sin(theta);
    y_theta_r = -x*sin(theta)+y*cos(theta);
    
    for m = 0:S-1
        c2 = a^(-m);
        x_theta = c2*x_theta_r;
        y_theta = c2*y_theta_r;
        g = c0*exp(-0.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)+c1*x_theta+psi);
        G{m+1,n} = g;        %g_normalize(g);
    end
end
end %gabor_fn2

function g = g_normalize(g)
    % NORMALIZATION of positive and negative values to ensure that the integral of the kernel is 0.
    % This is needed when phi is different from pi/2.
    ppos = find(g > 0); %pointer list to indices of elements of result which are positive
    pneg = find(g < 0); %pointer list to indices of elements of result which are negative
    
    pos = sum(g(ppos));  % sum of the positive elements of result
    neg = -sum(g(pneg)); % abs value of sum of the negative elements of result
    meansum = (pos+neg)/2;
    
    if (meansum > 0)
        pos = pos / meansum; % normalization coefficient for negative values of result
        neg = neg / meansum; % normalization coefficient for psoitive values of result
    end
    
    g(pneg) = pos*g(pneg);
    g(ppos) = neg*g(ppos);
end